c     
c      Based on John Peacock's Software on Halo Model
c      Further modifying by T. Miyaji, Jul 2009 to separate the case where 
c      one gal is at the center. Also one can choose between 
c      NFW and Moore99 profile
c
c      Modified by T. Miyaji Jul 2008
c      For Cross-correlation function
c      
c      BUG Fix: 04-MAY-2009: in powint (rk) 
c            rm_crit_0 = min(rm_crit(1),rm_crit(2))  (max->min)
c
      real*8 x_pdat(1000),y_pdat(1000,2)
      common/cospar/om_m,om_l,om_b,h,p_index,sig8
      common/cosmo_today/om_m0,om_l0,sig80
      common/wtpar/rm_crit(2),alpha(2),isam !N(M) parameter for two samples
      common/look/rnu_look(1000),rvir_look(1000),con_look(1000)
      common/norm/pfac
      common/pdat/x_pdat,y_pdat,n_pdat
      common/xint/ih,rad 
      external f_wprpint, f_xiint
      character*60 comd
      character*40 argc1,argc2,cch1,cch2,rbin
      character*5 typ, bias_only
      real lgme1,lgme1s,lgme1c,lgme2,lgme2s,lgme2c
c
c     Reading Commad Line Argument
c
c     Arguments 1:z 2:N(M) model/read by separate module 3-4.[Cosmo]
c               5: bin 6: m99 or nfw, 7:bias_only
c
      call getarg (0,comd)
      narg=iargc()
      if (narg.ge.1) then
        call getarg (1,argc1)
      end if

      if (narg.lt.2.or.index(argc1,'-h ').eq.1) then
         call print_help()
         go to 999              ! end
      end if
c
      read (argc1,*) z 
      call getarg (2,argc2)
c     Arg 2 is specific to the N(M) models
c      read (argc2,*) alpha(1),rm_crit(1),alpha(2),rm_crit(2)
c     
c      if (rm_crit(1).lt .29) rm_crit(1)=10**rm_crit(1) !Log Input
c      if (rm_crit(2).lt .29) rm_crit(2)=10**rm_crit(2) !Log Input
c       
      if (narg.ge.3) then
         call getarg (3,cch1)
      else
         cch1 = '- '
      endif
c
      if (narg.ge.4) then
         call getarg (4,cch2)
      else
         cch2 = '- '
      endif        
c
      call init_csm(cch1,cch2)
      sig8=sig80/grow(z)
c
      if (narg.ge.5) then
         call getarg (5,rbin)
         if (index(rbin,'- ').le.0) read (rbin,*) rmin,rmax,nr
      else
         rbin = '- '
      endif        
      if (index(rbin,'- ').ge.1) then
         rmin = 0.05
         rmax = 80.
         nr = 40
      endif
      dlnr = (log(rmax)-log(rmin))/float(nr)
c
      if (narg.ge.6) call getarg (6,typ)
      if (typ.ne.'nfw'.and.typ.ne.'NFW') typ = 'M99'        
c
      bias_only = 'No'
      if (narg.ge.7) call getarg(7, bias_only)
c
      om_m1=omega_m(z)
      om_l1=omega_v(z)
      om_m=om_m1
      om_l=om_l1      
      pfac=1.
c
      pfac=sig8/sigint(8.)
      pfac=pfac*pfac
c
      call init_nm (argc2)  ! Initialize Model calc defined in a separate module.
      call init_look(z) ! Initialize Lookup table for Power Spectrum?
c
c      write(0,'(A)') '! Initilaizations done'
c      Get Header Info on Bias Paramater and 
c
      b_lin1s = bias(1,0)
      b_lin1c = bias(1,1)
      b_lin1 = bias(1,-1)
      b_lin2s = bias(2,0)
      b_lin2c = bias(2,1)
      b_lin2 = bias(2,-1)
c     Effective Mean Mass and model galaxy density
      call meff (lgme1s,gden1s,1,0)
      call meff (lgme1c,gden1c,1,1)
      call meff (lgme1,gden1,1,-1)
      call meff (lgme2s,gden2s,2,0)
      call meff (lgme2c,gden2c,2,1)
      call meff (lgme2,gden2,2,-1)
c
c     Create Power spectrum 1 and 2 halo terms
c     
      write (6,'(A,A)') '! DMH Profile: ',typ
      write (6,'(A)') '! Samp.   bias  b_sat  b_cen LgMeff'//
c                        12345 123456 123456 123456 123456
     $     '  sat    cen    Density      sat       cen  '
c           123456 123456 123456789 123456789 123456789 
      write (6,'(A,6(1X,F6.2),3(1PE10.2))') '! Samp1',b_lin1,b_lin1s,
     $     b_lin1c,lgme1,lgme1s,lgme1c,gden1,gden1s,gden1c
      write (6,'(A,6(1X,F6.2),3(1PE10.2))') '! Samp2',b_lin2,b_lin2s,
     $     b_lin2c,lgme2,lgme2s,lgme2c,gden2,gden2s,gden2c      
c 
      if (bias_only(1:1).eq.'Y'.or.bias_only(1:1).eq.'y'.or.
     $     bias_only(1:1).eq.'1') goto 999  ! Do not calc CFs, display
      ! bias, meff, dens only
c
      call mk_pdat(b_lin1*b_lin2,typ) ! param: bias^2.
      write(0,'(A)') '! Power spectrum table made'
c          
      write (6,'(A)') '!  r,r_p      w_p(r_p)[h^-1 Mpc]'//
     $ '               xi(r)'
      write (6,'(A)') '! [h^-1Mpc]   1-halo     2-halo      sum  '//
     $     '     1-halo     2-halo       sum   '
c                       1234567890  1234567890 1234567890 1234567890 
c               1234567890 1234567890 1234567890
      rat_1h2h=1.0
      wp_1h_p = -1.0
      xi_1h_p = -1.0
      wp_2h_p = -1.0
      xi_2h_p = -1.0
      do i=0,nr
         r=rmin*exp(i*dlnr)
         if (rat_1h2h.gt.1e-4) then
            ih=1    ! ih is in common block passed onto integrand/1-halo term
            wp_1h=max(xiint(f_wprpint,r),2e-29) 
!                    Integrate Powerspec to wp(rp)
            xi_1h=max(xiint(f_xiint,r),2e-29)   
!                     Integrate Powerspec to xi(r)
         else
            wp_1h=1e-29
            xi_1h=1e-29
         endif
         ih=2 ! 2-halo term
         wp_2h=xiint(f_wprpint,r)
         xi_2h=xiint(f_xiint,r)
         rat_1h2h = xi_1h/xi_2h
         write (6,'(F9.3,2X,6(1X1PE10.2))') r,wp_1h,wp_2h,wp_1h+wp_2h,
     $        xi_1h,xi_2h,xi_1h+xi_2h
      enddo 
      write (6,'(A)') '! Created by '//comd(1:index(comd,' '))//
     $     argc1(1:index(argc1,' '))//argc2(1:index(argc2,' '))//
     $     cch1(1:index(cch1,' '))//cch2(1:index(cch2,' '))//
     $     rbin(1:index(rbin,' '))//typ
 999  end

      
c
      subroutine print_help()
      character*80 com_name
      call getarg (0,com_name)
      print '(A)',''
      print '(A)',' Usage:'//com_name(1:index(com_name,' '))//
     $     'z N(M)_pars [cosmo_par1 [cosmo_par2 [rbin [prof]]'
      print '(A)',''
      print '(A)','   N(M)_pars: model paramaters (see below)'
      print '(A)',''
      print '(A)','   cosmo_par1: h,OmeM,OmeL (default=0.7,0.3,0.7)'
      print '(A)','   cosmo_par2:'//
     $     ' sigma8_0,Ome_b,p (default=0.85,0.045,0.97)'
      print '(A)','   rbin:'//
     $     ' rmin,rmax,nbin (log binning/default=0.05,80,60)'
      print '(A)','   prof:'//
     $     ' m99 (default) or nfw'
      print '(A)',''
      print '(A)','Lengths are in h{-1} Mpc, masses are in h{-1} Msun.'
      print '(A)',''
      call model_help ()
      return
      end
c
